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The dynamic aeroelastic instability of flutter is an important factor in the design of modem 
high-speed, flexible aircraft. The current trend is toward the creative use of composites to delay 
flutter. To obtain an optimum design, we need a accurate as well as efficient model. As a first 
step towards this goal, flutter analysis is carried out for an unswept composite box beam using 
a linear structural model and Theodorsen’s unsteady aerodynamic theory. Structurally, the wing 
was modeled as a thin-walled box-beam of rectangular cross section. Theodorsen’s theory was used 
to get 2-D unsteady aerodynamic forces, which were integrated over the span. A free-vibration 
analysis is carried out using the theory of Ref. [4]. These fundamental modes are used to get the 
flutter solution using the V-g method. Future work is intended to build on this foundation. 


Structural Preliminaries 


A thin- walled box beam is considered for the analysis. A coordinate system x = {£ i x 2 £3} is 
defined along the undeformed wing with X\ along the elastic axis. 

The displacements, rotations, strains and curvatures are denoted by u, 6, 7 and k respectively 


U = i 

M 

u 2 

>6= < 


► 'y = < 

7u 

2712 ' 

► k = i 

2 


l“3. 


1^3 J 


l 2713 J 


1*3 J 


( 1 ) 


lYansverse shear deformations are negligible, so that 712, 713 = 0 - Thus 
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' o, ] 


6= < 

— u' 3 > and/c = < 

-«5 


i <4 J 1 

t4 


( 2 ) 

(3) 


Using the theory developed by Badir et.al. (1992) we get a constitutive equation of the form 
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where F is the force along the wing and Mi, M2, Mj are the torsional and two bending moments. 
The [C] 4x 4 is the cross sectional stiffness matrix and Is obtained U6ing the theory. 
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Energy Formulation and Raleigh Ritz Method 

Langrange equations of motion can be written in terms of the kinetic energy (T) and the strain 
energy (U) as 


dqi 


Kinetic energy is given by 
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where m is the mass per unit lenth, f is a matrix of center of mass offset from the elastic axis and 
i is the matrix of polar moment of inertia. For our problem Eq. (6) beoomes 
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Strain energy is given by 
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Now generalized coordinates are introduced to make it a finite dimensional system. We represent 
the deflections it and 6 with finite number of modeshapes yielding 
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Here n,o,p,g are the number of extentional, torsional, vertical bending, and inplane bending modes, 
respectively. Thus, we also have 
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Substituting Eq. (10) and Eq. (12) in the kinetic and strain energy expressions given by Eq. (7) 
and Eq. (9), the Lagrange’s equations can be expressed as 
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Calculation of {Q} will be discussed later. 

For the present analysis, n extensional, o torsional, p vertical bending and q inplane bending un- 
coupled modes are taken as {<£ U i}’s, {^,}’s, {dual’s and (</» ua }’s respectively. Analytic expressions 
for the modes were taken from Ref. [2]. Analytic expressions for the above integrals were derived 
for the case of uniform cross section beams. 


Free vibration analysis 

For free vibration analysis, {<5} = 0, thus Eq. (14) becomes 

wm + {f<m = 0 a?) 

Now assuminq simple harmonic motion, i.e. {9} = {gje****, we get an eigenvalue problem, whose 
eigenvalues give the frequencies and eigenvectors the modeshapes of vibration. 

[*]{«> = «*[«]{»} (18) 

Given below is a comparison of free-vibration frequencies of a Bending-Twist Coupled Beam with 
experimental values. A graphite epoxy cantilever plate beam is considered. The results validate 
the free-vibration analysis code. 
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Coordinate Transformations 

Aeroelastic analyses using simple structural models have traditionally relied upon aerodynamic 
variables such as angle of attack (a) and plunging motion (h) to describe the aerodynamic forces on 
the airfoil. Because the present analysis incorporates a more sophisticated structural model, it was 
found to be more convenient to express the aerodynamic parameters in terms of the wing kinematic 
variables. 

The wing section has an angle of attack a in a uniform freestream of velocity U. The reference 
frame d is oriented such that U = -Ua 2 (see Figure 1). The reference frame b is oriented with 
respect to the undeformed wing section and is centered at the section referenoe line. Thus, it is the 
same as the structural coordinate system x, i.e. = Xj. The referenoe frame B is a small angle 
transformation with respect to b caused by deformation of the section such that 

f Bi 1 [1 u 2 U3I f h ) 

Ib 2 \= -u' 2 1 0, { 62 \ (19) 
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The section of the undeformed beam is rotated at a steady-state angle of attack a 0 from the 
freestream such that the transformation from a to b is 

h] 1 0 ° I [ a, 1 

i b 2 > = 0 cosa 0 sina 0 j a 2 > (20) 
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In this analysis, we assume that the steady-state angle of attack is small (i.e. the same order of 
magnitude as the angles in Eq. 19). This allows us to make the approximations 

cos a 0 =l and sina 0 = a 0 . (21) 

Combining Eq. (19), Eq. (20) and Eq. (21), and discarding the higher order terms, the transfor- 
mation between the freestream reference frame and the deformed beam sectional reference frame 

r 1 w 2 u ' 3 ] rar 

< B 2 > = — u 2 1 + 0i < a 2 V (22) 

(Ba J L-U3 — (a 0 + 0i) 1 J la 3j 

This transformation shows the section angle of attack as a function of time to be 

a(t) = a 0 + 6i(t) (23) 

A 

The deflections Ui along the undeformed reference axes bi are also expressed in terms of the 
freestream reference frame a. Carrying through the small angles assumption of Eq. (21) and as- 
suming the displacements also to be small, it can be verified that 

UiXi = Uibi « Ujdj (24) 

This is a gross approximation, but it is allowed in this analysis because the small angle and 
small displacement assumptions are within the limitations of Theodorsen’s theory. From Eq. (24), 
the section plunging motion is 

h=-u 3 . (25) 
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Development of Generalized Aerodynamic Forces 

Lagrange’s equations were developed from kinetic and strain energy formulations for a uniform, 
cantilevered, thin-walled, closed-section beam and are given by 


Mq + Kq = Q 


(26) 


where M is the n x n generalized mass matrix, K is the n x n generalized stiffness matrix, Q is the 
n x 1 generalized foroe vector, and q are the generalized coordinates. There are two components 
of the generalized forces that must be considered in the flutter analysis: the aerodynamic loadings 
and the structural damping. This analysis uses 2-dimensional incompressible strip theory to obtain 
the aerodynamic loadings. The effect of structural damping is introduced artificially using the V-g 
method, as will be discussed in a later section. 

The vector of generalized aerodynamic forces is developed in from the principle of virtual work. 
The virtual work done on a two-dimensional airfoil section is expressed as: 


W = (L + D) • 8r + M • 6$ (27) 

where L, D and M are the section lift, drag and moment, respectively; <5r is the virtual displacement 
vector, and 8& is the virtual rotation vector: 
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The aerodynamic loadings are calculated with respect to the freestream reference frame a and 

— ♦ — #■ — V 

about the midchord, such that L = La 3, D = —Da 2 and M = Ma\. The reference line offset from 
the midchord is baB 2 . The virtual work can be rewritten using Eq. (28) in terms of the virtual 
displacements and the virtual rotation about the reference line: 


W= P t 8u + R t 6u' (29) 

where P and R are 4x1 column matrices and 
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The virtual work can be expressed in terms of generalized coordinates, using the expansion 


u = $q 


(30) 


where 4> is the modal matrix and q is the generalized coordinate vector. The virtual work in terms 
of generalized coordinates is 


W = P T ^6q + R T ~6q 
aq aq 


or, from Eq. (30), 


SW = (P r $ + R T V)8q 


5 



The coefficient of 6q is the transpose of the generalized aerodynamic forces on the section. Inte- 
grating over the wing span, we obtain the generalized aerodyanmic forces as 

Q = [ L ($ T P + V T R)dx (31) 

Jo 

The elements of P and R are, in terms of the lift, drag and moment, 

Pj =0 Pi =0 

P2 = M + baL 4- boD(ct 0 4- 9i) R% = 0 ^21 

P 3 = —L R 3 = -Mu ' 2 K ] 

P 4 = —D Pi = baDx 4 — baLu' 3 — Mu' 3 


2-Dimensional Strip Theory (Theodorsen) 

Expressions for the unsteady aerodynamic lift and moment on an airfoil were initially developed by 
Theodorsen. Assuming simple harmonic motion the total lift and moment can be expressed as 

L = Lq 4- Le^ 1 /oo\ 

M - Mo + Me" 1 y ) 


here Lq and M 0 are the steady components of the lift and moment associated with the steady- 
state angle of attack a 0 ; L and M are the complex magnitudes of the unsteady lift and moment, 
respectively; and u is the natural frequency. The steady components of the aerodynamic loadings 
are given from Ref. [1] in terms of reduced frequency k = ub/U as 
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In these equations, p is the air density and b is the semichord. Because this flutter analysis is but 
a first attempt using a primitive aerodynamic theory, the steady lift and moment coefficients are 
assumed to be those of a thin symmetric airfoil, i.e. 


Ci 0 = 27r q 0 and Cm 0 = na 0 


(35) 


The unsteady components of the lift and moment are derived from the complex magnitudes of 
the lift and moment coefficients in Ref. [1] for the case of simple harmonic pitching and plunging 
motions. These unsteady loadings are: 

L = -npb 3 u 2 {Lh^ + [L a - (^ + a)Pfc]a} 

M = 7tp6V{[M A - j + \M„ - - (i + a)(M h - jt k )]a} (36) 

The coefficients L a , L*, M a and Mh are functions of reduced frequency k and the Theodorsen 

function C(k), and are given from Ref. [1] as 

L a = 5 — r(l + 2C(fc)) — jpjC(fc) M q = f-£ 

L h = 1-f C{k) M h = | ^‘ } 
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The total lift and moment can be written in terms of the kinematic variables by recognizing 
that, from Eq. (23) and Eq. (25) 

h = — U 3 a — di 

Also, in accord with the assumption of simple harmonic motion, 

it3 = Use** 61 — Oxe** 


Expressions for the total lift and moment are thus 

L _ 27r ^ _— a 0 - 7rpi 3 a; 2 { - + [L a - (- 4- a)L/,]0i} 

M = ^a 0 + ^V{[M k -iL,]^5 (38) 
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The drag on the airfoil section is developed by introducing the small perturbations in velocity 
and angle of attack into the equation 

D = p6f/ 2 C£>(a) (39) 

From the geometery and noting the approximations made in Eq. (24) and Eq. (23), the velocity is 
U + u 2 and the angle of attack is o = a 0 + 0i- The velocity U is expressed in terms of reduced 
frequency k. In addition, it is recognized that, for simple harmonic motion, 

• * _ it. it • 
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Substituting these perturbations into Eq. (39) and using a Taylor series expansion, the linear 
approximation of the unsteady section drag is: 
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The generalized unsteady aerdynamic forces are obtained by substituting Eq. (38) and Eq. (40) 
into the relations for P and R in Eq. (32), and subsequently integrating Eq. (31). Neglecting 
higher-order terms, the non-zero elements of P and R become 
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The natural frequency term u> 2 can be factored out of each term in the above expressions, and these 
forces can be rewritten in the matrix form 


P = u 2 Au + u 2 B R = u ?Du' 


(42) 


The 4x4 matrices A and D and the 4x1 column vector B are shown in the appendix. Using 
the modal expansion in equation (30) and substituting the above expressions into Eq. (31), the 
generalized aerodynamic forces become: 

Q = u> 2 j ($ T A$ + & t D&) dx| q + u 2 £ $ t B dx (43) 


Equations of Motion 

Analysis of the flutter problem is greatly simplified when the mass matrix M in Eq. (26) is a 
diagonal matrix. The mass matrix in our analysis, however, contains off-diagonal terms that account 
for inertial couplings in the beam that arise when the center of mass does not coincide with the 
reference line. To solve this problem, the equations of motion in Eq. (26) can be recast in terms of 
a diagonalized mass matrix. To do this we consider the problem of the wing in free vibration 

Mq + Kq = 0 (44) 

The eigenvectors ft and eigenvalues u> 2 of this problem are obtained and used to simplify the 
equations of motion in the flutter problem. The mass and stiffness matrides can be pre- and 
post-multiplied by the n x n eigenvector matrix T , where 

T=[?i g 2 ••• Qn ] 


such that 

T t MY q + Y T KYq = 0 

The matrix M D = T r A7T is the diagonalized mass matrix, and the matrix K D = Y T KY is the 
diagonalized stiffness matrix. From the solution of equation (44), it can be seen that 

K d = a > 2 r QM d 


Here we have introduced a reference frequency ujr that will ultimately simplify the flutter equa- 
tion. This reference frequency is arbitrary, and is traditionally taken to be the natural frequency 
of the wing in the first torsion mode. The n x n matrix Q contains the free vibration natural 
frequencies in the following form: 


n = 
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This diagonalization technique is applied to the equations of motion for the flutter problem. 
Using the eigenvectors and eigenvalues from the free vibration problem of Eq. (44), Lagrange’s 
equations (26) become 

M°q + u%SlM D q - u 2 Eq = u 2 t T £ $ T B dx (45) 

where E is the generalized aerodynamic forces from Eq. (43) premultiplied by the transpose of the 
eigenvector matrix such that 

S = T T [ L (4> r A$> + 4>TD4>) dx 
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Structural Damping 

As pointed out in Bisplinghoff, Ashley and Halfman [1], a flutter analysis cannot provide accurate 
results without accounting for the effects of structural damping in some fashion. In flutter specifi- 
cally, structural damping is the only mechanism through which energy is removed from the wing. A 
precise description of the structural damping in the wing is difficult to obtain; however, since it is 
small compared to the aerodynamic loadings, it is possible to approximate the structural damping 
by an artificial damping force. To characterize this artificial damping force, consider a single degree 
of freedom system undergoing the simple harmonic motion 

— _ -u»4 

x = x 0 e 


A viscous damping force opposing this motion could be 

F d = -/^ = -fiuxoe 1 " 1 
at 

The viscous damping force in this example lags the position vector by 90° and has an amplitude 
that depends on the frequency. The structural damping in the flutter problem is known to be 
frequency independent, so the viscous damping model is not a valid approximation in toto. It does 
show, however, the presence of a phase shift between the motion and the corresponding damping 
force. 

In the motion of a beam element in free vibration, an artificial structural damping can be 
associated with each elastic restoring force and torque by multiplying the restoring force by a small, 
artificial damping coefficient g { associated with each motion. In terms of generalized coordinates, 
this artificial damping force can be defined from Eq. (44) as 
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The coefficients g { have absolutely no physical meaning. The flutter condition occurs when one 
of these coefficients goes to zero in the flutter equation. Because the structural damping effects are 
small in comparison with the aerodynamic forces, we can for the purpose of this analysis assume 

that 

gi ~ g? ~ ~ g n cz g 

The individual values of g will be recovered from the solution of the flutter determinant, as will 
be seen. 


Flutter Equation 

The generalized artificial damping force vector can now be included in Eq. (45) to obtain the 
complete equations of motion for a wing in flutter. 

M°q + ^oj 2 R QM D q + u 2 R QM D q - u 2 E q = u 2 T T £ * T B dx (47) 

This is an inhomogeneous ordinary differential equation in q. The inhomogenous term can 
be dropped if we consider the generalized displacement to be composed of steady-state and time- 
dependent components, such that, assuming again simple harmonic motion, 

q = q 0 + qe iu>t 


Dividing the remaining terms by the natural frequency u 
flutter equation becomes 

a M D (ij|) (1 + ig) - M D - E 1 


and rearranging, 
9 = 0 
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The the flutter conditions are found by setting the flutter determinant to zero for a range of 
reduced frequencies k R and plotting the natural frequency u and artificial damping coeffient g with 
respect to the reduced velocity, l/k R . Flutter occurs at the reduced velocity where g goes to zero. 


The flutter determinant is 


det 



(1 + ig) - M 


d _ c 


(49) 


The components of the flutter determinant can be arranged in a convenient form that allows the 
determinant to be solved in terms of a single complex polynomial. This is done by introducing the 
parameter d 2 r , which is a given natural frequency of vibration for a specific mode (usually taken to 
be the first torsion mode). The components of this determinant are: 


[( s )’ 2 - 0 + #)] lari=j 

-Eij for i j 


(50) 


where 


The natural frequency u and the artificial damping g are recovered from Z for each given reduced 
frequency. 
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Preliminary Results and Future Research 

A flutter analysis was carried out for a Uniform Cantilever Beam. The result is compared with that 
obtained by Goland Ref. [3]. The flutter velocity obtained by the code written was 305 mph while 
Ref. [3] obtained 385mph. The discrepancy in the result may be attributed to some differences in the 
aerodynamic model used and some slight numerical errors in Ref. [3]. The formulation developed 
is quite general and can form basis for further work in this field. The structural model could be 
extended to include nonlinearities. The aerodynamic model could be updated to a more accurate 
doublet lattice method. To do a time-domain analysis for accurate results away from the flutter 
point, a finite-state aerodynamics could be utilized. 


Appendix 

The elements of the matrices A, B and D contained in equation (42) are presented below. 
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